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In a recent series of papers, Gebremariam, Bogner, and Duguet derived a microscopically based 
nuclear energy density functional by applying the Density Matrix Expansion (DME) to the Hartree- 
Fock energy obtained from chiral effective field theory (EFT) two- and three-nucleon interactions. 
Due to the structure of the chiral interactions, each coupling in the DME functional is given as 
the sum of a coupling constant arising from zero-range contact interactions and a coupling function 
of the density arising from the finite-range pion exchanges. Since the contact contributions have 
essentially the same structure as those entering empirical Skyrme functionals, a microscopically 
guided Skyrme phenomenology has been suggested in which the contact terms in the DME functional 
are released for optimization to finite-density observables to capture short-range correlation energy 
contributions from beyond Hartree-Fock. The present paper is the first attempt to assess the ability 
of the newly suggested DME functional, which has a much richer set of density dependencies than 
traditional Skyrme functionals, to generate sensible and stable results for nuclear applications. The 
results of the first proof-of-principle calculations are given, and numerous practical issues related 
to the implementation of the new functional in existing Skyrme codes are discussed. Using a 
restricted singular value decomposition (SVD) optimization procedure, it is found that the new 
DME functional gives numerically stable results and exhibits a small but systematic reduction of 
our test x 2 function compared to standard Skyrme functionals, thus justifying its suitability for 
future global optimizations and large-scale calculations. 

PACS numbers: 21.10.-k,21.30.+y,21.60.Jz 



I. INTRODUCTION 

One of the fundamental challenges of nuclear theory is 
to predict properties of nuclei starting from the underly- 
ing vacuum two- and three-nucleon interactions. While 
impressive progress has been made in extending the lim- 
its of ab-initio methods beyond the lightest nuclei [ll-Q , 
the nuclear energy density functional (EDF) approach 
is the only computationally feasible many-body method 
capable of describing nuclei across the mass table [J]. 
Driven by interest in the coming generation of radioac- 
tive isotope beam facilities, along with studies of astro- 
physical systems such as neutron stars and supernovae 
that require controlled extrapolations of nuclear prop- 
erties in isospin, density, and temperature, there is a 
large effort currently under way to develop nuclear en- 



* stoitsovmv@ornl.gov 
t kortclainenc@ornl.gov 
t bogncr@nscl.msu.edu 
§ thomas.duguct@cca.fr 
^ furnstahl.l@osu.edu 
** gcbrcmar@nscl.msu.edu 
TT |schu ncknf@or nl.gov| 



ergy functionals with substantially reduced global er- 
rors and improved predictive power away from stabil- 
ity. The Universal Nuclear Energy Density Functional 
(UNEDF) SciDAC-2 collaboration is one such effort that 
aims to develop a comprehensive theory of nuclear struc- 
ture and reactions utilizing the most advanced computa- 
tional resources and algorithms available, including high- 
performance computing techniques to scale to petaflop 
platforms and beyond [5[ . 

Well-known empirical Skyrme and Gogny EDFs are 
typically characterized by 10-15 coupling constants ad- 
justed to selected experimental data. Despite their sim- 
plicity, such functionals provide a remarkably good de- 
scription of a broad range of bulk properties such as 
ground-state masses, separation energies, etc., and to a 
lesser extent of certain spectroscopic features of known 
nuclei. They are also widely employed, with some sucess, 
in studies of complex nuclear phenomena such as, e.g., 
large-amplitude collective motion. However, their phe- 
nomenological nature often leads to parameterization- 
dependent predictions and does not offer a clear path 
toward systematic improvements. 

One possible strategy is to provide microscopic con- 
straints on the analytical form of the functional and the 
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values of its couplings from many-body perturbation the- 
ory (MBPT) starting from the underlying two- (NN) 
and three- nucleon (NNN) interactions [a-fl2l|. Recent 
progress in evolving chiral effective field theory (EFT) 
interactions to lower momentum using renormalization 
group (RG) methods [HO] (see also [H, d|) is ex- 
pected to play a significant role in this effort, as the 
Hartrec-Fock approximation becomes a reasonable (if not 
quantitative) starting point. This suggests that the the- 
oretical developments and phenomenological successes of 
EDF methods for Coulomb systems may be applicable to 
the nuclear case for low-momentum interactions. 

However, even with these simplifications, the MBPT 
energy expressions are written in terms of density ma- 
trices and propagators folded with finite-range interac- 
tion vertices, and are therefore non-local in both space 
and time. In order to make such functionals numeri- 
cally tractable in heavy open-shell nuclei, it is necessary 
to develop simplified approximations, for example based 
on the use of local densities and currents. At lowest or- 
der in MBPT (i.e., Hartree-Fock), the den sity matrix ex- 
pansion (DME) of Negele and Vautherin [20| provides a 
convenient framework to approximate the spatially non- 
local Fock energy as a local Skyrme-like functional with 
density-dependent couplings. This novel density depen- 
dence of the couplings is a consequence of the finite- 
range of the vacuum interactions, and is controlled by 
the longest-ranged components. Consequently, the DME 
can be used to map physics associated with long-range 
one- and two-pion exchange interactions into a local EDF 
form that can be implemented at minimal cost in exist- 
ing Skyrme codes. The rich spin and isospin structure 
of such interactions should improve the quantitative pre- 
dictive power of EDFs while in the same time retaining 
the connection of these functionals with the underlying 
microscopic theory of nuclear forces. 

We do not expect dramatic changes for bulk nuclear 
properties due to the tendency of pions to average out 
in spin and isospin sums, but we do expect interesting 
consequences for single-particle properties (which phe- 
nomenology tells us are sensitive probes of the tensor 
force) and systematics along long isotopic chains (which 
should be sensitive to the isovector physics coming from 
pion-exchange interactions). Another potentially signif- 
icant advantage of the DME functional is that two very 
different microscopic origins of spin-orbit properties (i.e., 
short-range NN and long-range NNN spin-orbit interac- 
tions) are treated on an equal footing. This is in con- 
trast to empirical Skyrme and Gogny functionals, where 
the density-independent spin-orbit coupling is consistent 
with the short-range NN spin-orbit interaction, but has 
no obvious connection with the sub-leading (but quanti- 
tatively significant) long-range NNN sources of spin-orbit 
physics. Although it is beyond the scope of the present 
paper, let us mention that a clear priority of future stud- 
ies will therefore be to examine if the DME-bascd func- 
tional is able to improve two major shortcomings of stan- 
dard Skyrme phenomenology [2l| : (i) the destructive in- 



terplay between tensor and spin-orbit terms that com- 
promise spin-orbit splittings and the evolution of nuclear 
shell with isospin, and (ii) the too-high location of high- 
l centroids compared to \ow-l ones, which compromise 
the shell position even for nuclei near the stability valley. 
Any positive change regarding these two points would sig- 
nificantly impact the performance and predictive power 
of EDF calculations dedicated to spectroscopy. 

Recently, Gebremariam et al. have used an improved 
formulation of the DME jUJ to construct a non-empirical 
Hartree-Fock energy functional from un-evolved chiral 
EFT two- and three-nucleon interactions through next- 
to-ncxt-to-leading-order (N 2 LO) [H, [HJ. The struc- 
ture of the EFT interactions implies that each cou- 
pling in the DME HF functional can be written as the 
sum of a density-independent (Skyrme-like) coupling con- 
stant arising from the zero-range contact interactions 
and a density-dependent coupling function arising from 
the long-range pion-exchange interactions. As discussed 
in Section [Hi in the present approach the separation 
of long- and short-distance physics at the HF level is 
used to motivate a semi-phenomenological functional. 
In particular, the DME coupling functions arising from 
finite-range pion-exchanges are not modified, while the 
density-independent couplings associated with the con- 
tact interactions are released for optimization to infinite 
nuclear matter and finite nuclei properties in order to 
mimic higher-order short-range correlation energy con- 
tributions. 

It is expected that the semi-phcnomcnological DME- 
bascd functional should perform at least as well as em- 
pirical Skyrme functionals because one still fits the same 
Skyrme coupling constants to data, the only difference 
being that the new EDF contains additional parameter- 
free coupling functions derived from the finite-range NN 
and NNN interactions. However, due to the highly 
non-trivial density-dependence carried by the DME cou- 
plings, there is no a priori guarantee that the imple- 
mentation will not be plagued with numerical instabil- 
ities or other technical difficulties that invalidate the 
approach. Consequently, the main goal of the present 
paper is to perform "proof-of-principle" calculations in 
which we (i) give the practitioner's view of how the DME 
functional can be implemented in existing Skyrme codes, 
and (ii) perform a restricted SVD optimization ("pre- 
optimization" ) of the density-independent couplings to 
verify that the new microscopically guided phenomenol- 
ogy does no worse than standard Skyrme functionals, 
thus justifying its suitability for future global optimiza- 
tions and large-scale EDF calculations. 

In our initial investigation, we constrain zero-range vol- 
ume parameters of the DME functional to reasonable 
values for equilibrium characteristics of infinite nuclear 
matter (INM), while zero-range surface parameters are 
obtained from a restricted optimization procedure using 
SVD techniques based on 72 even-even nuclei binding 
energies, and 8 odd-even mass (OEM) differences (4 neu- 
tron and 4 proton) . Our analysis will show that while the 
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DME-bascd functional is indeed more susceptible to in- 
stabilities, it can still be made sufficiently stable to carry 
out an optimization procedure for nuclei throughout the 
nuclear mass chart. It should be stressed that a detailed 
comparison of the quality of the DME functional against 
standard Skyrme predictions, and a fortiori experimen- 
tal data, is premature before applying a more rigorous 
global optimization. 

The rest of the paper is organized as follows. In Sec- 
tion nn we review how the DME can be used to map 
ab-initio MBPT energy expressions into the form of a 
local EDF and motivate the semi-phenomenological ap- 
proach used in the present work. The explicit form of 
the functional is given in Section IIII1 and the free pa- 
rameters entering the volume part of the functional are 
expressed in terms of infinite nuclear matter (INM) equi- 
librium characteristics in Section lTVl The restricted SVD 
optimization procedure used to fix free surface parame- 
ters in the functional is described in Section [Vj and a 
comparison of selected nuclear properties calculated with 
both the DME functional and the standard Skyrme func- 
tional is made in Section IVI1 Conclusions are given in 
Section EU while various formulas and technical details 
are collected in the Appendices. 



II. MICROSCOPICALLY MOTIVATED 
FUNCTIONAL 

A. DME exchange energy functional 

At lowest order, the Fock energy computed from un- 
evolved chiral interactions exhibits spatial non-localities 
due to the convolution of finite-range form factors with 
non-local density matrices. The central idea of the DME 
is to factorize the non-locality of the one-body density 
matrix (OBDM) by expanding it into a finite sum of 
terms that are separable in relative r = r\ — r 2 and 
center of mass R = (ri + r2)/2 coordinates. Adop ting 
notations similar to those introduced in Refs. [HjIIJ, one 
expands the spin-scalar and spin-vector parts (in both 
isospin channels) of the density matrix as 

™max 

Pt(ri,r 2 )^J2u„(kr)V n (R), (1) 

m max 

s t (ri,r 2 ) « n m (fcr) Q rn {R) , (2) 

where k is an arbitrary momentum that sets the scale for 
the decay in the off-diagonal direction, whereas II„(fcr) 
denotes the so-called II— functions that depend on the 
particular formulation of the DME, see Refs. [U l23| . 
In the present work, k is chosen to be the local Fermi 
momentum related to the isoscalar density through 

/o 2 \V3 

k = k F (R) = (— p (R)\ , (3) 



although other choices are possible that include ad- 
ditional r- and Ap-dependencies [25j . The functions 
{V n (R), Q m {R)} denote various local densities and their 
gradients { Pt (R),T t (R), J t (R),V p t (R), Ap t (R)}, which 
for time-reversal invariant systems are defined by 

p t (R) = pt(ri,r 2 )\ ri =r 2 =R , (4) 
T t (R) = Vi • V 2 pt(ri,r 2 )\ ri =r 2 =R , (5) 

J t (R) = -l(Vi - V 2 ) x S t (r 1 ,r 2 )\ ri =r 2 =R , (6) 

where the isospin index t = {0, 1} labels isoscalar 
and isovector densities, respectively. For example, the 
isoscalar local density is the sum po = p n + p p of neutron 
p n and proton p p densities, while the isovector local den- 
sity is the difference p\ = p n — p P - Analogous expressions 
hold for T t , Jt, and all other quantities labelled by the 
index t = {0, 1}. 

Applying the expansion in Eqs. H}-© to the non- 
local exchange (Fock) energy gives a spatial integral over 
a sum of bilinear and trilinear products of local densities. 
For time-reversal invariant systems (truncating the ex- 
pansion to second-order in gradients), the NN exchange 
energy becomes 

^ NN W« E fdn{g pp p 2 + g pT Ptr t + g pAp p t Ap t 
t=o,i ^ 

+ gi Vp J t .V Pt +g-l J J^, (7) 

while the NNN contribution yields 

E™[p] » J dR^gPopl + g^p pj + g^plro 

+ g p ' To pU + g p ° p ^ Po pin + g p ° Ap °p 2 Ap 
+ g^pjApo + gP^ A P-p p 1 Ap 1 

+ g paJ °pvJl + g paJ? PoJf + 9 PlJoJl PiJo ■ Ji 
+ .9 poVpoJ VoVp - Jo 

+ g po * plJl PoVpi • Ji + g Pl ^ poJl piVp • Ji 
+ g Pl vpl Jo pi Vpi • Jo + g pl ^ Jo P lV ■ Jo 

+ g p2 ^- h p\V ■ Jo + g pap ^' h PoPiV ■ J^ (8) 

where, for simplicity, the R-dependcncc of the lo- 
cal densities and DME couplings has been omitted. 
The R-dependence (or equivalently, isoscalar density- 
dependence via Eq. of the couplings arises from the 
integral of the finite-range NN and NNN interactions over 
various products of the II-functions, e.g., 

ff f (R) ~ J drr 2 U p (k F r) U p (k F r) T?(r) , (9) 

where in this example T^* (r) is the central component of 
the exchange force V(r)P a P T . 

If the objective is to derive a fully microscopic and 
quantitative EDF free from any fitting to data, then the 
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purely non-empirical HF functional of Ref. |23j is inad- 
equate since un-evolved chiral interactions generate too- 
strong coupling between low and high momenta for HF to 
be a reasonable zeroth-order approximation. Moreover, 
it is known that even if the interactions are softened by 
evolving to low momentum, it is still necessary to go to at 
least 2nd-order MBPT to obtain a reasonable description 
of bulk properties of infinite matter as well as binding en- 
ergies and charge radii of closed-shell nuclei. 

Unfortunately, a consistent extension of the DME pro- 
cedure beyond the Hartree-Fock level of MBPT has not 
yet been formulated. At this point in time, any attempt 
to microscopically construct a quantitative Skyrme-like 
EDF must therefore inevitably resort to either some ad 
hoc approximations (e.g., neglecting state-dependent en- 
ergy denominators) when applying the DME to iterated 
contributions beyond the HF level, and/or to the re- 
introduction of some phenomenological parameters to be 
adjusted to data [Til l20l. l26l - [28l ] . An example of the latter 
approach has been recently proposed in Refs. |22l |23| . 

B. Semi-phenomenological DME functional 

Schematically, the EFT NN and NNN potentials have 
the following structure 

v EFT = v« + Kt , (io) 

where V v denotes finite-range pion-exchangc interactions 
and V c t denotes scale-dependent zero-range contact terms 
encoding the effects of integrated-out degrees of free- 
dom (e.g., heavier meson exchanges, high- momentum 
two-nucleon states, etc.) on low-energy physics. Con- 
sequently, each DME coupling in Eqs. 0^© decom- 
poses into a density-independent coupling constant aris- 
ing from the zero-range contact interactions, and a 
density-dependent coupling function arising from long- 
range pion exchanges, e.g., 

= <?n*w+cr(Kt), (ii) 

and so on. Note that the zero-range V c t generates 
Hartree-Fock contributions that are identical in form to 
the standard Skyrme functional, which is hardly surpris- 
ing since such functionals were originally derived as the 
Hartree-Fock energy density resulting from a zero-range 
Skyrme "force" or pseudo-potential. 

Based on the clean separation between long- and 
short-distance physics at the HF level, a semi- 
phenomenological approach where the long-distance cou- 
plings (g™(R; Vk)) are kept as is, and the zero-range C t m 
are optimized to finite nuclei and infinite nuclear mat- 
ter properties was suggested in Refs. [lH, While this 
is an admittedly empirical procedure, it is motivated by 
the observation that the dominant bulk correlations in 
nuclei and nuclear matter are primarily short-ranged in 
nature, as evidenced by Brueckner-Hartree-Fock (BHF) 
calculations where the Brucckner G-matrix "heals" to the 
free-space interaction at sufficiently large distances. 



Therefore, while a Hartree-Fock calculation using un- 
evolved chiral NN and NNN interactions would provide 
a very poor description of nuclei, the application of the 
DME to such contributions nevertheless captures at least 
some of the non-trivial density dependencies that would 
arise from the finite-range tail of an in-medium ver- 
tex (e.g., a G-matrix or a pcrturbativc approximation 
thereof) that sums ladder diagrams in a more sophis- 
ticated many-body treatment. In this sense, one can 
loosely interpret the refit of the Skyrme constants to data 
as approximating the short-distance part of the G-matrix 
with a zero-range expansion through second order in gra- 
dients. In the following section, we describe how free pa- 
rameters entering the volume part of the proposed func- 
tional can be fixed to equilibrium properties of infinite 
matter, while the remaining free surface parameters can 
be fixed to properties of nuclei using a restricted SVD 
optimization procedure. 

III. IMPLEMENTATION OF THE DME 
FUNCTIONAL 

A. Notations 

We write the proposed semi-phenomenological DME 
functional 

E[p] = J H(r)dr , (12) 
in the following form 

w 

Htt>{r) = U&ptpv + U&ptTv + u£j t • J v 

+ U£ p pt Apt, + U t p t VJ Pt V • J v , (14) 

where the notation reflects that our attention is restricted 
to the ground states of even-even nuclei in the present 
paper. Consequently, only terms built out of time-even 
densities are shown explicitly. The density-dependence 
of the couplings has been omitted for brevity. Note that 
the strange "off-diagonal" isospin structure in Eq. (fT4")l 
is a consequence of absorbing an extra factor of po or 
pi into the definition of the U$ couplings, which allows 
the trilinear 3N contributions in Eq. ([5]) to be written 
in terms of the more familiar bilinear products of local 
densities, e.g., 

9 Ato p\t - {g plTo Pi}piT = U% Pi to , (15) 

and so on. 

The U™, couplings (where m runs over the bilinears 
{ptpv , ptTv , Jt • Jv , pt^pv , ptV • JV}) have the following 
general structure 

U%, = {C?+g?{u) + Po h?{u))8 tiV 

+ Pl ^(«)(1-V) , (16) 
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where u = fci?(R)/m 7r . The functions gl n (u) are ob- 
tained by applying the DME to the Fock-energy contri- 
butions from the finite-range pion-exchange parts of the 
chiral EFT NN interaction through N 2 LO. The functions 
htf,(u) and h™(u) originate from the finite-range part 
of the leading chiral NNN interaction (which appears at 
N 2 LO), and are related to the couplings in Eq. (jHJ by 

h? = Pog r° m , m £ {p 2 t ,p t r t ,p t Ap t ,...} , (17) 
hf t , = P ig pim 7 m e {p tPt ,,p t T t ,,p t Apt',...} .(18) 

The C t m parameters correspond to the zero-range V c t con- 
tributions, which as discussed in the previous section will 
be released for optimization. In this way, the proposed 
DME functional splits into two terms, 

E[p] = E ct [p]+E x [p] , (19) 

where the first term -E c t[p] = £?[p] g =/i=o collects all con- 
tributions from the contact part of the interaction plus 
higher-order short-range contributions encoded through 
the optimization to nuclei and nuclear matter, while 
the second term E^[p] = E\fJ[c n ,=a collects the long- 
range NN and NNN pion exchange contributions at the 
Hartree-Fock level. 

This leads to the following explicit form of the DME- 
based energy density 

U(r) = ^—t + H (r) + Hi(r) + H 2 (r) , (20) 
where, for t = {0, 1}, 

Ht(r) = (C£ + G( D pl + gf (u) + p hf {u))p 2 t 
+ {Cf+g^{u)+p Q h^{u))p t r t 
+ (C t p Ap + g pAp (u) + Po h pAp (u))p t Ap t 
+ (C pVJ + g pWJ (u) + p»h? J (u))p t VJ t 
+ {cf +gf{u) + p hf{u))j 2 (21) 

and 

% 2 (r) = PiKo( u )Pi r o + Pi h io P ( u )Pi^Po 

+ pihil{u) Ji J + pih^ J {u) Pl VJo ■ (22) 

The explicit forms of the functions <?™( M ) an d ^tv ( u ) nave 
been given in Ref.[23| and the companion Mathematica 
notebooks. In order to gain a feeling about the new den- 
sity dependencies entering through such couplings, we 
provide stripped-down "skeleton expressions" along with 
several explicit examples in Appendix [X] 

In order to facilitate the use of the DME functional 
in nuclear EDF calculations, we have also developed a 
general module written in FORTRAN 90 [Hj], which can 
easily be ported to any existing EDF solver. It contains 
all of the lengthy expressions for the DME couplings C/ t ™ , 
Eq. (fT6|) , their functional derivatives with respect to the 
density matrix, and numerically stable approximate ex- 
pressions at small u. The module also has the capability 
to calculate related infinite nuclear matter properties. 



B. Contact part 

The contact part -EctI/?] has the form of the standard 
Skyrmc functional 

n ct (r) = ^T + H$(r)+H?(r) ) (23) 

where 

Uf(r) = (C& + c( D pl)p\ + Cr Pt r t + C pAp Pt A Pt 

+ C pVJ p t \7J t + C t j2 J?. (24) 

This is illustrated in Appendix [B] where the link be- 
tween the coupling constants and the historical (t„,x„)- 
paramctcrization of Skyrme "forces" is explicitly given. 

As for the standard Skyrme functional, Eq. (f2"3")l con- 
tains 13 parameters, 

{c£,cf D ,C? Ap ,Cr,Cf,CF J ,l} , (25) 

which are to be released for optimization to infinite mat- 
ter and finite nuclei properties. While these parameters 
have exactly the same form as in the standard Skyrme 
functional, the existence of the long-range part in the 
functional will obviously modify their optimized values. 

C. The parameter 7 

Early versions of the Skyrme functional motivated the 
p 1 term appearing in Eq. (|24[) as arising from a zero- 
range NNN force, in which case 7 = 1. However, this 
interpretation was soon found to be problematic, as 7 = 
1 yields too large an incompressibility Q. Subsequent 
Skyrmc paramcterizations largely cured this difficulty by 
letting 7 float, with values typically between 1 /6 and 1/3. 

In EFT studies of dilute Fermi systems interact- 
ing with zero-range interactions, one finds similar non- 
integer powers of p appearing in the energy density, which 
can be traced to correlation (i.e., beyond HF) effects j30j . 
Even in this much simpler model system, where a con- 
trolled and well-defined EFT expansion is possible, it is 
interesting to note that one finds multiple non-integer 
powers of p occurring at low orders in the expansion. 
Given that the nuclear many-body problem is much more 
complicated, with many additional possible sources of 
non-analytic behavior due to the interplay of finite-range 
NN and NNN interactions and short-range correlation ef- 
fects analogous to those found in the dilute fcrmion sys- 
tem, the single non-integer p 1 term in Eq. (|23[) is proba- 
bly not justified on microscopic grounds. 

Nevertheless, in the short term we follow standard 
practice with a single p 1 term in the functional, which 
in our case can be loosely viewed as parameterizing the 
HF contribution of the NNN contact term, plus higher- 
order correlation effects that are implicitly included in 
the refit to data. However, ultimately one would like to 
revisit this issue to see if MBPT can be used to provide 
insight regarding the form of such non-analytic terms. 
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D. Finite-range part 



The finite range part E^lp] follows from 

U„{r)=HZ{r)+Ul{r)+UZ{r) , 
where, for t = {0, 1}, 



(26) 



(gf(u)+p hf(u))p 2 

+ {g pT ( u ) + po^ T {u))p t T t 



[gr p (u) + p h pAp (u))p t Ap t 
(gf(u)+ Po hf(u))j? 



+ (g pVJ (u) + p h pVJ (u))p t WJ t , (27) 



oVJ 



and 



+ Pl h(o(u)JiJ + Pl h p ^ J {u)p 1 VJ . (28) 

Couplings entering -E^p] are entirely determined in 
terms of the finite-range NN and NNN interaction pa- 
rameters, and are therefore frozen during the optimiza- 
tion procedure. In the present work, the values for the 
couplings that enter the finite-range chiral EFT interac- 
tions arc taken from Ref. 13111. 



E. Hartree NN contributions 

In general, it is possible to apply the DME to both 
Hartree and Fock energies so that the complete Hartree- 
Fock energy is mapped into a local functional. It is 
known since the original work of Negele and Vautherin, 
however, that treating the Hartree contributions exactly 
provides a better reproduction of the density fluctuations 
and the energy produced from an exact HF calculation 
[26l . [32( 1 . Restricting the DME to the exchange contribu- 
tion significantly reduces the self-consistent propagation 
of errors (26|. Moreover, treating the Hartree contribu- 
tion exactly generates no additional complexity in the 
numerical solutions of the resulting self-consistent equa- 
tions compared to applying the DME to both Hartree and 
Fock terms. Lending further support to Negele and Vau- 
therin's conclusions, we find that the present DME-bascd 
functional becomes extremely susceptible to numerical 
instabilities when the DME approximation is used for 
the Hartree terms, which immediately disappear when 
the Hartree terms are treated exactly. 

In the present work, a simplification is used such that 
finite-range NN Hartree contributions are treated in the 
Local Density Approximation (LDA). An exact treat- 
ment of these contributions and their influence on the 
results will be examined in a future investigation of the 
DME-based functional. 



IV. CONSTRAINING THE VOLUME TERM 
PARAMETERS 

A. INM with the DME-based functional 

We turn now to a discussion of how equilibrium prop- 
erties of infinite nuclear matter (INM) can be used to fix 

the 7 free volume parameters {C t p , C^, C t pT , 7} in the 
functional. In INM, the total energy per particle defines 
the saturation curve W(p n , p p ). Its derivation discards 
the Coulomb energy and all gradient terms that are zero 
for a homogeneous system, and substitutes the kinetic 
energy density with its Thomas-Fermi expression, which 
is exact in this case. Assuming spin saturation, one also 
disregards the terms involving the spin-orbit current den- 
sity, J t . 

The expansion of W(p n ,p P ) around the equilibrium 
density p c in a Taylor series in p = p n + p p and / = 
{Pn~ Pp)/P yields 

W(p, I) = W{p) + S 2 (p)I 2 + S 4 (p)I 4 , (29) 

pNM pNM t-NM 

W( P ) = — + —S P + — (dp) 2 , (30) 



A 

c I \ NM . 

£>2( P ) = a sym + 



p; 

£NM 

3p c 



*P+^(6pf 



18 P 2 



(31) 



-NM _NM 



where Sp = {p-p c ), while E /A, P NM , K 
L NM and AK NM denote the total energy per nucleon 
at equilibrium, the pressure, the nuclear matter incom- 
pressibility, the symmetry energy coefficient, the density 
derivative of af^, and the isovector correction to the in- 
comprcssibility at saturation density p c of nuclear mat- 
ter, respectively. The quartic term in / can be safely 
neglected in Eq. (|29| in practice. 

The INM equation of state (EOS) following from the 
functional Eq. (fT4|) takes the form 

W(I,P) = ^T 

+ {csl + c( D p 1 + gf («) + p K 2 («))p 

+ {Cf + cf D p^ + g( (u) + p h( (u))l 2 p 
+ + g p T (u) + pK T {u) + I 2 pKl{u))r 

(32) 



(Cf +g p l T (u) + ph pT {u))lT 1 



where 



37L! 
2 



1/3 



(33) 



r = ic P 2 / 3 ((l + /) 5 / 3 + (l-/) 5 / 3 ) , (34) 
r 1 = ic P 2 / 3 ((l + /) 5 / 3 -(l-/) 5 /3) , (35) 



5 V 2 



2/:! 



(36) 



Our strategy is to express unknown volume parameters 
in terms of nuclear matter equilibrium quantities p c , 
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E NM /A, K NM , a™, L NM , m*, and to*,, where m* and 
to* are the isoscalar and isovector effective masses, re- 
spectively. This strategy has been recently applied in the 
context of the optimization of pure Skyrme functionals 



B. Symmetric nuclear matter constraints 

Symmetric nuclear matter (SNM) is characterized by 
equal neutron and proton densities p n = p p = p/2, where 
p is the isoscalar density. All isovector terms are thus 
zero. The isoscalar kinetic energy density per particle 
from Eq. (03} is 



T = C P 



2/3 



(37) 



to* through 

v _ 27(7 + - 9L NM + 20(2 - 3 7 )C pr Pc r c 



C 



10 



CP' 



27 7 p c 

£((97 - 6)to* 1 - 12 7 + 5)r c + A 10 (t 



+ 



27 7 p c 



27a™ + 9L NM - 40C p > c t c 



7+1 



+ 
fi 2 



27jp 



7+1 



(42) 



(43) 



T = ^K"* - KT 1 ) 1 - 0i» - K»Pe (44) 
Ira p c 

where Cq T has already been determined by Eqs. 



2 2 



Parameters Cq , Cq^, and Cg T are then expressed in 
terms of E NM /A, p c , and m* through 



(3-(2-3 7 )to:- 1 )t c 
(38) 



2 1 F NM fi 2 

=^-[^ + 1 )V-|- 

3 7 p c A 2m 

+ A 00 (u c )] , 

2 1 pNM fc2 

+ A 0D (u c )] , (39) 
<T = U- K" 1 - 1) 1 - .9o PT K) - ^ r K)p c , (40) 



where r c and u c are the kinetic energy density and the 
dimensionless Fermi momentum at the saturation den- 
sity pc, respectively. Explicit expressions for A t o(u) and 
A t D{u) are given in Appendix [Cl 

As for the parameter 7, one can either leave it as a 
free parameter or eliminate it using the incomprcssibility 
if NM . The resulting expression is 

= -A' NM - 9£ NM /A + £(4TOr 1 - 3)r c + A 7 (u c ) 
7 9^ NM /A + 3f r (2TO*- 1 -3)r c + J B 7 K) 

(41) 

where explicit expressions for A~ t (u c ) and i? 7 (u c ) are 
given in Appendix [Cl 



C. Asymmetric nuclear matter constraints 

In the regime of isospin asymmetric INM, neutron and 
proton densities are different (po = p, pi = Ip) and the 

isovector terms contribute. One can therefore express 

2 2 

parameters Cf , C P D and Cf r in terms of a^™, L NM and 



D. Reference SLy4 properties 

In this work, we will often use as benchmark the SLy4 
parametrization of the Skyrme force [H[ . The optimiza- 
tion protocol of this interaction included data obtained 
from ab initio calculations in nuclear matter of Ref . [35| , 
and we see in Fig. [1] that the saturation curves obtained 
with SLy4 agree well with the ab initio results. We there- 
fore take the INM equilibrium characteristics of the SLy4 
parameterization as an acceptable set of values to be used 
in Eqs. (j38j)- (1131) . 



E 



NM 



.4 



K 



NM 



-15.97 MeV, 
229.9 MeV, 



0.1595 far 



^ = 1.44/1.25 , (45) 



a™ = 32.0 MeV, 



L™ = 45.96 MeV. 



The resulting parameters C t p , C P D , C pT and 7 are com- 
pared with the original SLy4 parameters in Table HI and 
the associated INM curves are shown in Fig. Q] for sym- 
metric matter (top panel) and pure neutron matter (bot- 
tom panel). 

As seen from Table HI the values of the refitted param- 
eters differ substantially from the original SLy4 values. 
The original value of 7 = 1/6, for example, increases to 
about 7=1/3 when only NN contributions are taken into 
account, but becomes almost equal to one when both NN 
and NNN contributions are accounted for at the N 2 LO 
level. Nevertheless, the EOS for the DME-based func- 
tional are practically identical to the original SLy4 curves 
for symmetric (top panel of Fig. [1} and pure neutron 
matter (bottom panel of Fig. [1} for densities relevant to 
nuclei. 

At higher densities (p > 0.3 fm~ 3 ) the symmetric mat- 
ter EOS (Fig. [T] top) remains completely predetermined 
by the imposed equilibrium values Eq. (|45p with a slight 
deviation towards the reference points when NNN con- 
tributions are taken into account (N 2 LO). Deviations in 
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FIG. 1. (Color online) INM saturation curves for symmet- 
ric (top panel) and neutron (bottom panel) nuclear matter 
calculated with the standard SLy4 Skyrme functional (dot- 
ted lines), and the DME functional at LO (dashed lines), 
NLO (dot-dashed lines), N 2 LO without NNN contributions 
(short dashed lines), and N 2 LO (solid lines). Ab-initio re- 
sults [35|] are plotted as reference points. All functionals are 
constrained to reproduce the reference INM saturation values 
from Eq. (f45|) . 



this density range become more visible for the neutron 
matter EOS (Fig. Q] bottom) , as isovector properties de- 
pend on the EFT order and/or whether the NNN force 
is taken into account (N 2 LO) or not (N 2 LO NN ). 

If one does not impose the nuclear incomprcssibility 
value A" NM = 229.9 MeV but instead varies the value 
of 7, one can trace the influence of the DME contribu- 
tions on the INM incomprcssibility A" NM . Such results 
are shown in Table [TT1 for two values of 7: 7=1/6, cor- 
responding to the SLy4 parameterization, and 7=1, 
originally proposed for Skyrme functionals |3(| ■ 

Contributions from the NN interaction generally re- 
duce the value of K NM by about 10-20 MeV. Conversely, 
NNN contributions give too-high values for K NM unless 
7 ~ 1, which brings the N 2 LO A NM value to the phys- 



TABLE I. Parameters entering the volume part of the stan- 
dard Skyrme functional and the DME functional at LO, NLO, 
N 2 LO without NNN contributions (N 2 LO NN ), and N 2 LO, 
all reproducing the same values Eq. (|45|) for INM saturation 
properties. 



Parameters 

IF 

2 

Wo 

2 

C 



SLy4 



LO 



NLO N 2 LO 



NN 



N 2 LO 



(~;p 



-933.342 

830.052 

861.062 

-1064.273 
57.129 
24.657 



■727.093 

477.931 

612.114 

■705.739 
33.885 
32.405 



-757.689 

477.931 

628.504 

-694.665 
18.471 
92.233 



0.16667 0.30622 0.287419 



-777.805 

677.296 

641.017 

-952.381 
26.0411 
-51.8352 
0.275049 



■607.108 

331.438 

■1082.85 

4383.27 
322.4 
-156.90 
1.06429 



TABLE II. Nuclear matter incompressibility iv NM (in MeV) 
calculated at two different values of 7. 



SLy4 



LO 



NLO N LO 



N 2 LO 



7 : 
7 : 



1/6 229.90 
1 356.36 



208.49 
336.34 



211.42 
338.96 



213.34 
340.68 



440.50 
244.98 



ically accepted range of 220-250 MeV. Interestingly, the 
fact that the preferred value of 7 is rather close to 1 
when finite-range NNN contributions arc explicitly ac- 
counted for seems to contradict the naive argument re- 
called in Sec. MI Cl that 7 should encode both the Hartree- 
Fock contribution of the NNN contact term (7=1) 
but also higher-order correlation effects producing non- 
integer values for 7. 



E. First test on surface parameters 

The equilibrium INM properties allow us to constrain 
the zero-range volume parameters of the DME-bascd 
functional (Table H} , but not the parameters C t pA ' 3 and 



a 



pVJ 



entering the surface part of the functional or the 



tensor parameters C t as their associated terms are zero 
in spin symmetric nuclear matter. We note right away 
that all tensor terms from both pion exchanges and the 
Skyrme-like contact terms are omitted in the present 
proof-of-principle investigation. The reasons for that are 
briefly discussed in Sec. IV Al 

As for a first test on the surface parameters of the 
DME-based functional, we keep the volume ones at the 
values estimated from INM (Table [J, and simply set 
Cj^ Ap and Cf VJ equal to zero. That is, we let the long- 
range DME part of the functional (£^-[p]) generate all of 
the surface contributions. Results of such calculations at 
LO for two benchmark nuclei 40 Ca and 208 Pb are shown 
in the second column of Table IIIII The comparison with 
the results by SLy4 (the first column in Tabic IIII[) shows 
unacceptable over-binding of about 140 MeV and 320 
MeV in 40 Ca and 208 Pb, respectively. 
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TABLE III. Comparison of calculated SLy4 energies (in MeV) 
with the results from DME calculations at LO for nuclei 40 Ca 
and 208 Pb: kinetic energy for neutrons T n and protons T p , 
volume energy Ev, surface energy Es (see also Table [IVJ and 
total energy Et-- Volume DME parameters are taken from 
Table HI1 Surface DME parameters C* t pAp and Cf J are set 
equal to zero (third column) or to their SLy4 values (fourth 
column). 







r pAp and n pVJ — n 


liUIll D-l-jV^i 




SLy4 


DME:LO 


DME:LO 






40 Ca 




T n 


321.788 


401.334 


313.497 


T P 


313.215 


393.286 


304.782 


E v 


-1161.116 


-1362.847 


-1130.494 


E s 


111.046 


7.889 


109.992 


Et 


-344.227 


-480.770 


-332.150 






208 pb 




T n 


2527.799 


2819.082 


2490.574 


T„ 


1336.341 


1481.145 


1321.254 


E v 


-6514.517 


-7092.892 


-6429.071 


E s 


315.116 


65.083 


316.181 


E T 


-1635.160 


-1950.690 


-1599.798 



As a second test case, we compute the LO DME re- 
sults in which we have taken the SLy4 values for Cf Ap 
and Cf VJ (third column on Table HIT]) . In this case, the 
results are much closer but now an under-binding is seen 
of about 10 MeV and 30 MeV in 40 Ca and 208 Pb, respec- 
tively. 

Table Mil suggests that it should be possible to opti- 
mize the surface parameters in the DME functional in 
a manner similar to the optimization done for standard 
Skyrme functionals. Broadly speaking, one could think 
of procedures based on semi-infinite NM properties, or 
on the leptodermous expansion of the functional. Both 
approaches would fix the parameters entering the surface 
part of the functional on a set of well-defined surface co- 
efficients. Alternatively, one can probe and constrain the 
surface parameters using properties of finite nuclei. This 
is the path taken in the present work and described in 
Section El 



of-mass coordinate [37| . This freedom introduces a factor 
(a 2 — a + 1/2) in front of the surface terms proportional to 
pAp, with a ranging between zero and one. In quantum 
chemistry studies of molecular exchange energies, taking 
a = 0, which corresponds to expanding p(ri,r2) asym- 
metrically about rg, was found to be the optimal choice 
by a large margin [37|. For nuclei, we find that the quan- 
tum chemistry choice leads to overly strong surface con- 
tributions, resulting in heavy under-binding, which can 
only be compensated by substantially reducing the asso- 
ciated surface constants C t pAp . Such a reduction, how- 
ever, leads to numerical instabilities of the HFB solution 
for practically all nuclei studied. 

In our test calculations, we are able to avoid such dif- 
ficulties by using the original Negele-Vautherin choice 
(a = 1/2), which corresponds to symmetrically expand- 
ing p(ri,r2) about the center-of-mass coordinate R = 
(i"i + Y2)/2. This choice minimizes the long-range pAp 
contributions, and leads to stable results for the nuclei 
studied in the present paper. 

Another source of instabilities has been found with re- 
spect to the isovector behavior of the DME-based func- 
tional. For example, increasing the value of the symme- 
try energy parameter a^J, generates a functional which 
has an instability that cannot be compensated by the 
contact part of the functional. As a matter of fact, the 
N 2 LO functional displays systematic instabilities in fi- 
nite nuclei calculations when it is defined with the val- 
ues of Table HU However, by slightly modifying the 
values of INM characteristics to a™. = 30 MeV and 

ay ill 

i NM = 40 MeV (compare with values in Eq. 05])), the 
N 2 LO functional becomes stable enough to carry out the 
SVD optimization of the surface parameters as discussed 
in the next section. The resulting EOS in both symmet- 
ric and pure neutron matter are shown on Fig. [5] where 
one sees that the latter is much better behaved at high 
density as one increases the EFT order. 

In the calculations presented below, we have managed 
thus far to avoid such instabilities in our zcroth-order 
optimization ( "prc-optimization" ) to fix the volume and 
surface parameters of the DME-based functional. How- 
ever, a more sophisticated global optimization of the 
DME-based functional will ultimately require a more pre- 
cise and detailed analysis of its stability properties to rule 
out the most common problems [Hj]. 



F. Stability of the DME-based functional 

Our preliminary analysis of calculations using different 
sets of C™ parameters has shown that the DME func- 
tional is somewhat more sensitive to instabilities than 
standard Skyrme functionals, especially when the N 2 LO 
NNN contributions are taken into account. Nevertheless, 
we have been able to sidestep these issues thus far with 
minimal modifications. 

The DME procedure, for example, contains a freedom 
in the choice of the coordinate system when expanding 
the one-body density matrix with respect to the center- 



V. PRE-OPTIMIZATION 

A global optimization of a given EDF parameteriza- 
tion becomes a very involved procedure as soon as one 
extends the data set to include obscrvablcs beyond nu- 
clear binding energies. Such optimization procedures are 
expensive, as they require a large number of functional 
evaluations. It is always useful for the global optimiza- 
tion to have preliminary estimates for the values of the 
functional parameters. 

In the present section, we describe a particularly con- 
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FIG. 2. (Color online) INM saturation curves for symmetric 
(top panel) and neutron (bottom panel) nuclear matter cal- 
culated with the modified SLy4' Skyrme functional (dotted 
lines) and the DME functional at LO (dashed lines), NLO 
(dot-dashed lines), N 2 LO (solid lines) using the parameters 
from Table [TV] 



venient method to generate approximate estimates for 
the surface parameters constrained to binding energies 
and odd-even mass (OEM) differences. The method uses 
the explicit linear dependence of the functional on the 
parameters C™. The procedure can be seen as a partic- 
ular implementation of the optimization algorithm based 
on the regression analysis of Refs. [39l Hoi]. 



function is then defined as 



1 



X 



n c i - Tlx 



E 

i=i 



(46) 



where rid denotes the total number of experimental data 
points, n x the number of free parameters, E° xp the exper- 
imental values while Ei are the corresponding calculated 
ones. The weights a, are 2 MeV for binding energies and 
50 keV for pairing gaps. The same weights were used in 
Ref. 

The binding energy E A computed from a Skyrmc-likc 
EDF depends linearly on the coupling constants such 
that one can write 



Ea 



N 



Cke A 



(47) 



where the index k i— > {m, t] runs over all N combina- 
tions of m — {p 2 , pr, pAp, pVJ, J 2 } and isospin t = 0,1. 
Furthermore, e A collects all contributions to the energy 
that do not explicitly depend on the functional parame- 
ters {Cfc}, i.e., kinetic and Coulomb energies as well as 
contributions coming from finite-range pion exchanges. 
The terms e k A entering Eq. (|47|) are integrals multiplying 
the associated parameter Ck in the EDF. For example, 
the term associated with C/T is 



„Pa T o _ 



= / Po{r)T (r)dr 



(48) 



where the densities are computed for nucleus A. Strictly 
speaking, e\ depends on the value of the coupling con- 
stants. However, in close vicinity of the local minimum 
the behavior of the functional can be approximated to be 
linear in the couplings. 

Starting with a given initial set of coupling constants 
{C^}, one can perform EDF calculations for all the nu- 
clei belonging to the data set. This provides values of the 
integrals e A for every k and A. If some of the coupling 
constants are left out from the optimization process, i.e. 
n x < N, their contribution to the total energy can be 
transferred to the constant e A . Substituting coefficients 
e A in Eq. (|47[) and requesting the minimum of the x 2 
which now contains linearized energies, one ends up with 
a system of rid linear equations for the n x unknown pa- 
rameters {Ck}- 



A. SVD optimization procedure 

We first define the \ 2 function to be minimized. Due 
to the limitations of the optimization method, see below, 
the experimental data set is reduced to nuclear binding 
energies (i.e. masses) and pairing gaps. Except for proton 
radii, the data set is the same as in Ref. [33[. The \ 2 



Tlx 

£ C k e k Ai = E e *v -e%, (i = 1, . . . ,n d > n x ) . (49) 

k=l 

Solving the system of equations P9"j) with respect to {Ck} 
using the SVD method gives the solution for the mini- 
mum of the x 2 [] under the linear hypothesis. Due to the 
previously mentioned (small) non-linearities in e A , the 
process has to be repeated iteratively until the true local 
minimum is found. 
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The procedure described above can be easily extended 
to functionals containing a pairing contribution. In the 
present work we include a mixed delta pairing [HI, 



^ = i/(i-£)M + vD* 



(50) 



where p c = 0.16 fm~ 3 , and where V q , p q denote the pair- 
ing strength and pairing density for neutrons (q = n) 
and protons (q = p), respectively. In EDF calculations 
with pairing, one can calculate average pairing gaps for 
neutrons (A„) and protons (A p ) using 



1 - 



2po 



Pq( r )Pq( r ) > (51) 



where N n , and N p are the number of neutrons and pro- 
tons, respectively. 

It is straightforward to add the pairing parameters 
V„ and V p to the SVD optimization procedure above 
because Eqs. ([50)1 and (f5"Tj) are linear with respect to 
them. The pairing energy Eq. ([50)) should be added to 
Eq. (gTJ) , while Eq. (gSJ| should be extended with the re- 
quirement that the calculated (neutron or proton) pairing 
gaps Eq. (|5Tj) for selected nuclei should approximate the 
odd-even mass (OEM) difference defining the experimen- 
tal pairing gap I\ q 



a Ajexp _ a A, 

q q 



(j = l,...,M g ) 



(52) 



Besides the adjustment to binding energies and OEM 
differences, it is necessary to impose boundaries to the 
domain of variation of the free parameters due to the 
approximate nature of the nuclear functional and the in- 
complete set of experimental obscrvablcs used. For exam- 
ple, one cannot release for optimization parameters Cf^ J 
controlling spin-orbit contributions simultaneously with 
parameters C/ 2 governing tensor contributions [2l|. In 
the present proof-of-principle investigation, we drop all 
tensor contributions from pion exchanges entering the 
DME-bascd functional. That is, the amplitudes TJ^J are 
set to zero so that we can avoid optimizing the corre- 
sponding contact terms. Eventually, such terms will be 
subject to detailed investigation for the reasons explained 
in the Introduction. Even then, all remaining parame- 
ters C t m cannot be fully released for a SVD optimization. 
Such attempts lead to a negative isovector effective mass 
(experimental data cannot constrain this quantity [HI) 
or to a too-small INM saturation density (charge radii 
are not included in the present optimization). 



B. Preliminary parameterization 

The SVD optimization procedure has been performed 
with respect to the six (N = 6) parameters Cq Ap , Cf Ap , 
Cg VJ , C pVJ , V n , and V p using the binding energies of 30 
spherical and 42 deformed nuclei (i.e., M = 72), neutron 



pairing gaps of 4 nuclei (M n — 4), and proton pairing 
gaps of another 4 nuclei (M p = 4). Nuclear proper- 
ties have been calculated as in Ref. [33| using the HF- 
BTHO solver [43| in the Lipkin-Nogami regime of ap- 
proximate particle number projection [44| . For compari- 
son, the same SVD optimization procedure has also been 
performed for the standard Skyrmc functional starting 
from the SLy4 parameterization. The results define a 
new parameterization referred to as SLy4'. 



TABLE IV. Parameters entering the volume, surface, and 
pairing parts of the standard Skyrme functional and the 
DME-based functional at LO, NLO, and N 2 LO. The asso- 
ciated x 2 values and root-mean-square deviations (RMSD) 
are also shown. 



Parameter SLy4 SLy4' 



LO 



NLO 



N 2 LO 



Volume Parameters 

-933.342 -727.093 -757.689 -607.108 

830.052 474.871 477.931 316.939 

861.062 612.104 628.504 -1082.854 

-1064.273 -705.739 -694.665 -4369.425 

57.129 33.885 18.471 322.4 

24.657 32.405 92.233 -156.901 

0.16667 0.30622 0.287419 1.06429 



c 2 
c 2 

2 

r<P T 
C{ T 

7 







Surface Parameters 




°0 


-76.287 


-76.180 


-67.437 


-63.996 


-197.132 




15.951 


24.823 


21.551 


-9.276 


-12.503 




-92.250 


-92.959 


-95.451 


-95.463 


-193.188 


C pVJ 


-30.75 


-82.356 


-65.906 


-60.800 


37.790 



Pairing Parameters 
-258.992 -232.135 -241.203 -241.484 -272.164 
-258.992 -244.050 -252.818 -252.222 -286.965 



Vo 



SVD Optimization Results 

12.5002 2.1235 1.837 1.7662 1.7884 

7.008 2.6931 2.5539 2.5143 2.590 

0.1297 0.0828 0.0587 0.0554 0.0476 

0.094 0.0988 0.0902 0.0866 0.0706 



X 

(-B)rmSD 
(A n /RMSD 
(A p )rmSD 



The values of the parameters resulting from the SVD 
optimization, together with the \ 2 values and the result- 
ing root mean square deviations (RMSD), are shown in 
Table ITVl The first column in Table HVl shows the results 
for SLy4. Since our HFBTHO calculations are performed 
under the Lipkin-Nogami procedure, SLy4 leads to quite 
a high value of x 2 ps 12.5. The second column in Ta- 
ble QV] shows the results for SLy4'. The resulting x 2 is 
about six times smaller, \ 2 ~ 2.12. 

The optimal DME-based functional further reduces 
the value of x 2 , but by a much less significant amount. 
This is nevertheless a remarkable result keeping in mind 
the extremely involved structure, i.e. the rich density- 
dependence of the coupling functions, of the DME con- 
tributions and the fact that they do not contain optimiza- 
tion parameters. Interestingly, the N 2 LO parameteriza- 
tion (last column of Table IIV[) performs as well as stan- 
dard Skyrmc functionals with a reasonable incomprcss- 
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ibility K NM = 230 McV, even though its contact part is 
characterized by a density-dependent power 7 ~ 1. 

Because the DME-bascd functional is found to provide 
as good (or slightly better) a description of bulk prop- 
erties of finite nuclei as standard Skyrme functionals, a 
refined (and much more costly) global optimization can 
be undertaken that will eventually lead to a systematic 
analysis of its ability to improve the known deficiencies 
of standard functionals. Such a task is left for a future 
work. With a much less ambitious objective in mind, 
the following section gives some insight into the novelties 
that could be expected from DME functionals in typical 
nuclear structure applications by collecting a sample of 
results obtained with the parametrization of Table IIVI 



VI. SELECTED RESULTS IN FINITE NUCLEI 

This section summarizes results for selected nuclei, 
comparing their properties calculated with the standard 
Skyrme functional (SLy4' set of parameters) and the 
DME-bascd functional with the parameters given in Ta- 
ble IIVI Let us recall again that a detailed comparison of 
the DME functional to experimental data is premature 
until a more rigorous global optimization is performed. 

As the first example, Fig. [3| compares experimen- 
tal single-particle energies of [45j for the nucleus 208 Pb 
with the canonical single-particle energies calculated with 
SLy4' and the LO, NLO, and N 2 LO DME functionals. 
In general, the comparison shows that the DME func- 
tional docs not modify the Skyrme results significantly. 
The largest deviations are seen in the N 2 LO case mainly 
due to the stronger spin-orbit contact contribution (com- 
pare the values of Cf 7,1 from Tabic EJ). In |H] it was 
found that Skyrme functionals perform poorly on single- 
particle energies. The marked differences between SLy4' 
and N 2 LO single-particle spectra indicate that this situ- 
ation may be improved. 

However, as discussed in the previous section, all ten- 
sor contributions (contact term and finite-range contri- 
butions) have been set to zero to avoid the unnaturally 
large and strongly cancelling values for C t pVJ and Cf 
that arise in the present optimization protocol. Since the 
interplay between tensor and spin-orbit terms is crucial 
to understanding the evolution of nuclear shell structure, 
a detailed comparison of single-particle energies to data 
and standard Skyrme functionals is not appropriate with 
the present restricted optimization. 

Similar comparisons between SLy4', LO, NLO, and 
N 2 LO results are shown in Fig.2]for the two-neutron sep- 
aration energies (top panel), neutron rms radii (middle 
panel) and the average neutron gaps (bottom panel) for 
nuclei in the Ni chain in the region up to the neutron drip 
line. Again, the comparison in Fig. 0] shows similar be- 
haviors for the Skyrme and DME-based functionals, with 
the largest deviations coming at N 2 LO. Let us remem- 
ber, though, that small differences in separation energies 
can play an important role for reliable predictions of the 



position of the neutron drip line. As in the case of dif- 
ferent Skyrme parameterizations, the DME functionals 
could lead to a shift in this prediction. 

The fact that the Skyrme and DME functionals pro- 
duce very similar results for this pool of observables is 
in fact encouraging, since we do not want to lose the 
good features that phcnomcnological functionals based 
on the Skyrme force have acquired over the years. On 
the other hand, it would also be disappointing if the 
rich, microscopically-derived and non-trivial density de- 
pendence of the DME functionals could not bring in new 
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FIG. 3. (Color online) Comparison of 208 Pb experimental 
single-particle energies (in MeV) for neutrons and protons 
with canonical single-particle energies, calculated with the 
standard SLy4' functional and the LO, NLO, and N 2 LO DME 
functionals. 
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FIG. 4. (Color online) Two-neutron separation energies (top 
panel), neutron rms radii (middle panel), and average neu- 
tron pairing gaps (bottom panel) for the Ni isotopic chain of 
nuclei in the region of the neutron drip line, calculated with 
the standard SLy4' functional and the LO, NLO, and N 2 LO 
DME-based functionals. 



physics that can not be captured by Skyrme functionals. 

In this respect, nuclear deformation is an excellent 
probe, as it reflects the competition between the bulk 
properties of the interaction and its single-particle con- 
tent. Located right after the onset of deformation, the 
nucleus 100 Zr is characterized by the coexistence of 3 min- 
ima, oblate, spherical and prolate, the relative position of 
which is highly sensitive to the interaction. Fig. [5] shows 
that, in contrast to the Skyrme functional, the oblate and 
prolate minimum for DME functionals have almost the 
same energy (shape coexistence). At N 2 LO, the differ- 
ence is even more marked, as the spherical minimum dis- 
appears and is shifted at small prolate deformation. This 
behavior of the DME functionals can probably be related 
to a combination of small differences in single-particle en- 
ergies of closed-shell nuclei, see Fig. [31 as well as rather 
different surface bulk properties, see the value of cou- 
pling constants in Table IIVI Indeed, it is a particularity 
of that new generation of EDF parameterizations to pro- 
vide surface and spin-orbit terms with density-dependent 
couplings. 



FIG. 5. (Color online) Deformation energy for the nucleus 
100 Zr calculated with standard SLy4' functional (dotted line, 
squares) and the DME functional in LO (dashed line), NLO 
(dot-dashed line), and N 2 LO (solid line, circles). 



Another example of systematic differences is seen in 
the proton rms radii along the Ca isotopic chain as shown 
on Fig. |51 Since r.m.s. radii were not included in the 
pre-optimization of the DME functional, the particular 
value of the proton radius is irrelevant. However, the 
isotopic trend is a marker for the iso- vector channels of 
the functional, and the differences of slopes between the 
Skyrme and DME functionals, and curvature between 
LO, NLO and N 2 LO might be indicative of new physics. 




22 24 26 28 

Neutron Number N 

FIG. 6. (Color online) Comparison of the proton rms radii for 
nuclei along the Ca-chain calculated with the standard SLy4' 
functional (squares) and the DME functional in LO (dashed 
line), NLO (dot-dashed line), and N 2 LO (solid circles). 

Whether these changes improve or deteriorate the 
quality of the current functional with respect to the ex- 
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pcrimcntal data is irrelevant, as the paramctrization of 
Table [IV] should only be thought of as a prototype. A 
detailed study of the capabilities of DME functionals to 
reproduce experimental data, as well as more system- 
atic comparison with standard parametrizations of the 
Skyrme functional, will be performed once a more com- 
prehensive optimization procedure, such as the one used 
in Ref. [33[, has been carried out. 



VII. CONCLUSIONS 

In the present paper, we have given a practitioner's 
view of how the microscopically motivated DME func- 
tional of Gebremariam et al. (23j | . which possesses a 
much richer set of density dependencies than tradi- 
tional Skyrme functionals, can be implemented in exist- 
ing Skyrme HFB codes. Empirical infinite nuclear matter 
properties are used to constrain the volume parameters 
of the Skyrme-likc part, followed by a restricted singular 
value decomposition (SVD) optimization procedure to fix 
its density-independent surface parameters. We find that 
the proposed functional gives numerically stable results 
and exhibits a small but systematic reduction in x 2 com- 
pared to standard Skyrme functionals, thus justifying its 
suitability for future global optimizations and large-scale 
calculations. 

The DME-bascd functional takes the same general 
form as standard Skyrme functionals, with the key dif- 
ference that each coupling is composed of a density- 
independent Skyrme-likc piece that is optimized to data, 
plus a density-dependent coupling function determined 
solely (no free parameters) from the HF contributions 
of the underlying finite-range NN and NNN interactions. 
In this way, the functional is split into a parameter-free 
"long-range" piece that is directly linked to underlying 
NN and NNN pion-exchange contributions (treated at 
the HF level), and a short-range piece that is identical in 
form to the standard Skyrme functional with parameters 
that are optimized to data. 

After reviewing the structure of the DME-based func- 
tional and motivating the semi-phcnomenological ap- 
proach used in the present work, it was demonstrated 
how the free contact parameters entering the volume part 
of the functional can be eliminated in a one-to-one fash- 
ion in terms of equilibrium infinite nuclear matter char- 
acteristics. The influence of the finite-range DME con- 
tributions to symmetric and neutron INM was investi- 
gated with the most significant modification seen in the 
N 2 LO case when the three-body interaction is taken into 
account. In this case, we find a reasonable incompress- 
ibility K NM = 230 McV with a Skyrme parameter 7 ~ 1, 
a result which cannot be achieved within the standard 
Skyrme functional. 

A preliminary (pre-) optimization procedure for sur- 
face and pairing parameters using the SVD-optimization 
algorithm was performed using the binding energies of 
72 spherical and deformed nuclei, as well as 8 odd-even 



mass differences. It was found that the pre-optimized 
DME-based functional performs as well or slightly bet- 
ter than the standard Skyrme functional with respect 
to the optimized binding energies and OEM differences. 
The same is true also for other nuclear characteristics, 
e.g., nuclear rms radii, pairing gaps, separation energies, 
single-particle energies, as well as for nuclei that are not 
included in the optimization. These preliminary results 
are very encouraging, as they imply that more elabo- 
rate global optimizations of the DME functional will, at 
the very worst, preserve the already impressive level of 
performance provided by traditional Skyrme functionals, 
and very likely improve on them. 

The present results have been obtained under two re- 
strictions, which are not expected to modify our general 
conclusions and which will be removed in a future work. 
First, the NN Hartree contributions have been treated 
within the local density approximation. However, they 
can be easily taken into account exactly, as the compu- 
tational cost is the same as for the calculation of the 
Coulomb direct term, which is already included in nu- 
clear EDF calculations. 

Much more important, insofar as it plays a central role 
in future investigations of the spectroscopic and single- 
particle properties of the new functional, is the neglect 
of the tensor contributions in the present study. The 
issue is that in the present optimization procedure, the 
contact tensor coupling constant cannot be released for 
optimization together with the spin-orbit coupling con- 
stant, as the optimization drives both couplings to un- 
naturally large (and strongly canceling) values. It should 
be noted that this difficulty is not specific to the DME- 
bascd functional, as similar issues arise with the stan- 
dard Skyrme functional. The next step, removing both 
limitations mentioned above, is to apply a complete opti- 
mization procedure with the DME-based functional and 
perform systematic comparisons to the standard Skyrme 
functional and experimental data. Work in this direction 
is already in progress. 
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Appendix A: DME skeleton expressions 

The lengthy analytic expressions for the DME cou- 
plings tend to obscure their underlying structural sim- 
plicity. Therefore, it is more illuminating to display the 
couplings in a "skeleton form" that still conveys its non- 
trivial density dependence. 

The DME coupling g" l (u) is given as a sum of LO, 
NLO, and N 2 LO contributions (recall u = kp/m^), 
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where t = 0, 1 and the index m runs over the standard 
bilinear forms {p 2 , p t T t , ptApt, Jf , Jt^Pt}- These contri- 
butions are of the following generic form 
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where a 9 k 



+ 7l Vl + « 2 log(l + 2u 2 + 2uy/l + u 2 ) , (A3) 

5( u )In=lo = a 2+P2 lo g(! + « 2 ) + 72 arctan(u) , (A4) 

a fc( M )> K = PU U )> and 7fc = 7fe(«) are 
rational polynomials in it, with their dependence on t 
and m not explicitly shown. The explicit expressions 
for k = 0, 1,2 are given in Ref [23| and the companion 
Mathematica notebooks. 

In a similar way, the DME couplings h^, (it) that collect 
the N 2 LO NNN contributions read 

h£,(tt) = a h + $ log(l + 4u 2 ) + ft (log(l + 4u 2 )) 2 

+ 7q arctan(it) + 7^ (arctan(2it)) 2 



72 log(l + 4u 2 ) arctan(2w) , 



(A5) 



where the explicit expressions for the rational polynomi- 
als a\ = a%{u), p% = (3%(u), and 7^ = 7^ (it), with their 
dependence on t and m are not explicitly shown. 



Appendix B: Coupling constants and 
tx-parameterization 

Using the following explicit one-to-one relation be- 
tween C|" parameters and the (t,x) Skyrme parameters, 
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and substituting them into Eq. (j23p. one can verify that 
the contact part of the DME functional (|23|) is equivalent 
to the well-known Skyrme energy density: 
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where in neutron-proton notation q = (n,p), densities 
without an index stand for the total densities, e.g., p = 
Pn + p P , T = r„ + Tp, and J = J n + J p . 
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Appendix C: DME functions for INM 

The explicit expression for functions appearing in the 
INM equations are 
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where prime and double prime denote the first and second 
derivative with respect of u, respectively. 
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